files = ['HH';'HV';'VH';'VV'];
pol = ['RL';'RR';'LL';'LR'];
Ts = [11,20,30,40,50,55,60,65,70,80,100,130,160,190,220,250,275,296];
offsets=[]
%Ts = [80]%,100,130,160,190,220,250,275,296];
figure
hold on
EsVH=zeros(18,991);
dataVH=zeros(18,991);
for num=1:length(Ts)
    filename=['VH',num2str(Ts(num)),'.txt'];
    if Ts(num)==80
        filename=['VH',num2str(Ts(num)),'_2.txt'];
        data=readmatrix(filename);
        plot(data(26:end,1),data(26:end,2)+50*(num-1),'.-','Color',[(num)/20 0 0]);
        Es=data(26:end,1);
        Is=data(26:end,2);
        dataVH(num,:)=data(26:(26+990),2);
        EsVH(num,:)=data(26:(26+990),1);
    else
        data=readmatrix(filename);
        plot(data(260:end,1),data(260:end,2)+50*(num-1),'.-','Color',[(num)/20 0 0]);
        Es=data(260:1250,1);
        Is=data(260:1250,2);
        dataVH(num,:)=data(260:(260+990),2);
        EsVH(num,:)=data(260:(260+990),1);
        ylim([20,1000])
        title('VH')
    end
end
figure
imagesc(Es,Ts,dataVH)
colormap('jet')
clim([10,50])

files = ['HH';'HV';'VH';'VV'];
pol = ['RL';'RR';'LL';'LR'];
Ts = [11,20,30,40,50,55,60,65,70,80,100,130,160,190,220,250,275,296];
%Ts = [80]%,100,130,160,190,220,250,275,296];
figure
hold on
EsHV=zeros(18,991);
dataHV=zeros(18,991);
for num=1:length(Ts)
    filename=['HV',num2str(Ts(num)),'.txt'];
    if Ts(num)==80
        filename=['HV',num2str(Ts(num)),'_2.txt'];
        data=readmatrix(filename);
        plot(data(26:end,1),data(26:end,2)+50*(num-1),'.-','Color',[(num)/20 0 0]);
        Es=data(26:end,1);
        Is=data(26:end,2);
        dataHV(num,:)=data(26:(26+990),2);
        EsHV(num,:)=data(26:(26+990),1);
    else
        data=readmatrix(filename);
        plot(data(260:end,1),data(260:end,2)+50*(num-1),'.-','Color',[(num)/20 0 0]);
        Es=data(260:1250,1);
        Is=data(260:1250,2);
        dataHV(num,:)=data(260:(260+990),2);
        EsHV(num,:)=data(260:(260+990),1);
        ylim([20,1000])
        title('HV')
    end
end
figure
imagesc(Es,Ts,dataHV)
colormap('jet')
clim([10,50])
%% offset
%{
offsetsVH=zeros(18,1);
offsetsHV=zeros(18,1);
offseteqn='(a1*exp(-(x-b1)^2/c1^2)<64555)*a1*exp(-(x-b1)^2/c1^2)+(a1*exp(-(x-b1)^2/c1^2)>64555)*64555';
fo = fitoptions('Method','NonlinearLeastSquares',...
            'StartPoint',[0,0.2]);
ft = fittype(offseteqn,'options',fo);

for num=1:length(Ts)
    [curve2,gof2] = fit(EsVH(num,:)',dataVH(num,:)',ft);
    offsetsVH(num)=curve2.b1;
    [curve3,gof2] = fit(EsHV(num,:)',dataHV(num,:)',ft);
    offsetsHV(num)=curve3.b1;
end

for num=1:length(Ts)
    offsetsVH(num)=sum(EsVH(num,:).*(dataVH(num,:)>60000))/sum(dataVH(num,:)>60000);
    offsetsHV(num)=sum(EsHV(num,:).*(dataHV(num,:)>60000))/sum(dataHV(num,:)>60000);
end

offsetsHV-offsetsVH
%}
%% offset2
offsetsVH=zeros(18,1);
offsetsHV=zeros(18,1);
offseteqn='a1*exp(-(x-b1)^2/c1^2)+d1';
fo = fitoptions('Method','NonlinearLeastSquares',...
            'StartPoint',[10001,0,2,10]);
ft = fittype(offseteqn,'options',fo);


for num=1:length(Ts)
    EsVHfit2=EsVH(num,:);
    dataVHfit2=dataVH(num,:);
    EsHVfit2=EsHV(num,:);
    dataHVfit2=dataHV(num,:);
    EsHVfit2(dataHVfit2>10000)=300;
    dataHVfit2(dataHVfit2>10000)=0;
    EsVHfit2(dataVHfit2>10000)=300;
    dataVHfit2(dataVHfit2>10000)=0;
    [curve2,gof2] = fit(EsVHfit2',dataVHfit2',ft);
    offsetsVH(num)=curve2.b1;
    figure(22)
    plot(EsVHfit2',dataVHfit2','o')
    
    hold on
    plot(-10:0.01:10,curve2(-10:0.01:10))
    xlim([-10,10])
    hold off
    figure(23)
    [curve3,gof2] = fit(EsHVfit2',dataHVfit2',ft);
    plot(EsHVfit2',dataHVfit2','o')
    
    hold on
    plot(-10:0.01:10,curve3(-10:0.01:10))
    xlim([-10,10])
    hold off
    offsetsHV(num)=curve3.b1;
end

offsetsHV-offsetsVH

%% Gaussfit

EsHVfit=EsHV(:,870:960);
EsVHfit=EsVH(:,870:960);

dataHVfit=dataHV(:,870:960);
dataVHfit=dataVH(:,870:960);

eqn='a1*exp(-(x-b1)^2/c1^2)+d1';
fo = fitoptions('Method','NonlinearLeastSquares',...
            'StartPoint',[20,298,5,10]);
ft = fittype(eqn,'options',fo);
bsHV=[]
bsHVerr=[]
dsHV=[]
figure
hold on
for idx= 1:length(Ts)
    E=EsHVfit(idx,:);
    I=dataHVfit(idx,:);
    [curve2,gof2] = fit(E',I',ft);
    bsHV=[bsHV,curve2.b1];
    dsHV=[dsHV,curve2.d1];
    berr=confint(curve2,0.68);
    bsHVerr=[bsHVerr,berr(2,2)/2-berr(1,2)/2];
    plot(E,I+30*(idx-1),'.-')
    plot(285:0.1:310, curve2(285:0.1:310)+30*(idx-1))
end

figure(30);
errorbar(Ts,bsHV-offsetsHV',bsHVerr,'o')

bsVH=[]
bsVHerr=[]
dsVH=[]
figure
hold on
for idx= 1:length(Ts)
    E=EsVHfit(idx,:);
    I=dataVHfit(idx,:);
    [curve2,gof2] = fit(E',I',ft);
    bsVH=[bsVH,curve2.b1];
    dsVH=[dsVH,curve2.d1];
    berr=confint(curve2,0.68);
    bsVHerr=[bsVHerr,berr(2,2)/2-berr(1,2)/2];
    plot(E,I+30*(idx-1),'.-')
    plot(285:0.1:310, curve2(285:0.1:310)+30*(idx-1))
end
figure(30);
hold on
errorbar(Ts,bsVH-offsetsVH',bsVHerr,'o')
%% imagesc plot
figure
[X,Y] = meshgrid([EsHV(1,:)],[Ts,Ts]);
pcolor(X-[offsetsHV;offsetsVH],Y,[dataHV;dataVH])

shading flat
colormap jet
clim([10,50])
xlim([285,310])

figure
[X,Y] = meshgrid(EsHV(1,:),Ts);
pcolor(X-offsetsHV,Y,dataHV-dsHV')
shading interp
colormap jet
clim([0,40])
xlim([285,310])
ylim([0,300])


Tplot=0:305;
dataplot=zeros(length(Tplot),991*2);
offsetplotHV=[];
offsetplotVH=[];
Tline=[];
for num=1:(length(Ts)-1)
    Tline(num)=(Ts(num)/2+Ts(num+1)/2)-0.1;
end
for idx=1:length(Tplot)
    numT=sum(Tplot(idx)>Tline)+1
    dataplot(idx,:)=[dataHV(numT,:),dataVH(numT,:)];
    offsetplotHV=[offsetplotHV,offsetsHV(numT)];
    offsetplotVH=[offsetplotVH,offsetsVH(numT)];
end


figure
[X,Y] = meshgrid([EsHV(1,:),EsVH(1,:)],Tplot);
pcolor([X(:,1:991)-offsetplotHV',X(:,992:1982)-offsetplotVH'],Y,dataplot)

shading interp
colormap jet
xlim([285,310])
clim([10,50])
%% 
%{
figure(30)
hold on
bsall2=[bsVH-offsetsVH',bsHV-offsetsHV'];
Tsall2=[Ts,Ts];
a=2.0091;
c=300.2778;
a=1.939;
c=299.9;
eqnfit=@(x) -a*(1+2./(exp((c)./(2*0.69501*x))-1))+c;
plot(0:300,eqnfit(0:300))
a=2.2035;
c=300.7223;
a=2.117;
c=300.3;
eqnfit=@(x) -a*(1+2./(exp((c)./(2*0.69501*x))-1))+c;
bspart2=[bsVH(9:end)-offsetsVH(9:end)',bsHV(9:end)-offsetsHV(9:end)'];
Tspart2=[Ts(9:end),Ts(9:end)];
plot(0:300,eqnfit(0:300))
%ylim([290,300])

%eqnfit1='-a*(1+2/(exp((c)/(2*0.69501*x))-1))+c';

%}

%% seperate peaks
bs=[293,295,297,299,300.9,303]+0.5;
eqn='a1*exp(-(x-b1)^2/c1^2)+a2*exp(-(x-b2)^2/c1^2)+a3*exp(-(x-b3)^2/c1^2)+a4*exp(-(x-b4)^2/c1^2)+a5*exp(-(x-b5)^2/c1^2)+a6*exp(-(x-b6)^2/c1^2)+e*x+f';
sp=[10,30,50,30,10,10, bs(1), bs(2), bs(3), bs(4), bs(5), bs(6),.4,0,0];
bs=[293,295,297,299,300.9,303]+.54;
eqn2='a1*exp(-(x-b1)^2/c1^2)+a2*exp(-(x-b2)^2/c1^2)+a3*exp(-(x-b3)^2/c1^2)+a4*exp(-(x-b4)^2/c1^2)+a5*exp(-(x-b5)^2/c1^2)+e*x+f';
sp2=[10,30,50,30,10, bs(1), bs(2), bs(3), bs(4), bs(5),.4,0,0];
eqn3='a1*exp(-(x-b1)^2/c1^2)+a2*exp(-(x-b2)^2/c1^2)+a3*exp(-(x-b3)^2/c1^2)+e*x+f';
sp3=[10,20,20, bs(2)+.33, bs(3)+.3, bs(4)+.35,.4,0,0];


fo = fitoptions('Method','NonlinearLeastSquares',...
            'StartPoint',sp);
ft = fittype(eqn,'options',fo);
bsHV=[]
bsHVerr=[]
dsHV=[]
figure
hold on
for idx= 1:length(Ts)
    fo = fitoptions('Method','NonlinearLeastSquares',...
                'StartPoint',sp);
        ft = fittype(eqn,'options',fo);
    if Ts(idx)<45
        E=EsHVfit(idx,:);
        I=dataHVfit(idx,:);
        [curve2,gof2] = fit(E',I',ft);
        bsHV=[bsHV;[curve2.b1,curve2.b2,curve2.b3,curve2.b4,curve2.b5,curve2.b6]];
        dsHV=[dsHV,curve2.c1];
        berr=confint(curve2,0.68);
        bsHVerr=[bsHVerr;[berr(2,7)/2-berr(1,7)/2, berr(2,8)/2-berr(1,8)/2, berr(2,9)/2-berr(1,9)/2, berr(2,10)/2-berr(1,10)/2, berr(2,11)/2-berr(1,11)/2, berr(2,12)/2-berr(1,12)/2]];
        plot(E,I+30*(idx-1),'.-')
        plot(285:0.1:310, curve2(285:0.1:310)+30*(idx-1))
    end
    if Ts(idx)>45 && Ts(idx)<66
        fo = fitoptions('Method','NonlinearLeastSquares',...
                'StartPoint',sp2);
        ft = fittype(eqn2,'options',fo);
    
    E=EsHVfit(idx,:);
    I=dataHVfit(idx,:);
    [curve2,gof2] = fit(E',I',ft);
    bsHV=[bsHV;[curve2.b1,curve2.b2,curve2.b3,curve2.b4,curve2.b5,0]];
    dsHV=[dsHV,curve2.c1];
    berr=confint(curve2,0.68);
    bsHVerr=[bsHVerr;[berr(2,6)/2-berr(1,6)/2, berr(2,7)/2-berr(1,7)/2, berr(2,8)/2-berr(1,8)/2, berr(2,9)/2-berr(1,9)/2, berr(2,10)/2-berr(1,10)/2,0]];
    plot(E,I+30*(idx-1),'.-')
    plot(285:0.1:310, curve2(285:0.1:310)+30*(idx-1))
    end
    if Ts(idx)>66
        fo = fitoptions('Method','NonlinearLeastSquares',...
                'StartPoint',sp3);
        ft = fittype(eqn3,'options',fo);
        E=EsHVfit(idx,:);
    I=dataHVfit(idx,:);
    [curve2,gof2] = fit(E',I',ft);
    bsHV=[bsHV;[0,curve2.b1,curve2.b2,curve2.b3,0,0]];
    dsHV=[dsHV,curve2.c1];
    berr=confint(curve2,0.68);
    bsHVerr=[bsHVerr;[0, berr(2,4)/2-berr(1,4)/2, berr(2,5)/2-berr(1,5)/2, berr(2,6)/2-berr(1,6)/2,0,0]];
    plot(E,I+30*(idx-1),'.-')
    plot(285:0.1:310, curve2(285:0.1:310)+30*(idx-1))
    end
end

figure(40);
errorbar(Ts,bsHV-offsetsHV,bsHVerr,'o','Color',[0.85,0.33,0.10],'LineWidth',1,'MarkerFaceColor',[0.85,0.33,0.10],'MarkerSize',4)
ylim([290,310])


bsVH=[]
bsVHerr=[]
dsVH=[]
figure
hold on
for idx= 1:length(Ts)
    fo = fitoptions('Method','NonlinearLeastSquares',...
                'StartPoint',sp);
        ft = fittype(eqn,'options',fo);
    if Ts(idx)<45
        E=EsVHfit(idx,:);
        I=dataVHfit(idx,:);
        [curve2,gof2] = fit(E',I',ft);
        bsVH=[bsVH;[curve2.b1,curve2.b2,curve2.b3,curve2.b4,curve2.b5,curve2.b6]];
        dsVH=[dsVH,curve2.c1];
        berr=confint(curve2,0.68);
        bsVHerr=[bsVHerr;[berr(2,7)/2-berr(1,7)/2, berr(2,8)/2-berr(1,8)/2, berr(2,9)/2-berr(1,9)/2, berr(2,10)/2-berr(1,10)/2, berr(2,11)/2-berr(1,11)/2, berr(2,12)/2-berr(1,12)/2]];
        plot(E,I+30*(idx-1),'.-')
        plot(285:0.1:310, curve2(285:0.1:310)+30*(idx-1))
    end
    if Ts(idx)>45 && Ts(idx)<66
        fo = fitoptions('Method','NonlinearLeastSquares',...
                'StartPoint',sp2);
        ft = fittype(eqn2,'options',fo);
    
    E=EsVHfit(idx,:);
    I=dataVHfit(idx,:);
    [curve2,gof2] = fit(E',I',ft);
    bsVH=[bsVH;[curve2.b1,curve2.b2,curve2.b3,curve2.b4,curve2.b5,0]];
    dsVH=[dsVH,curve2.c1];
    berr=confint(curve2,0.68);
    bsVHerr=[bsVHerr;[berr(2,6)/2-berr(1,6)/2, berr(2,7)/2-berr(1,7)/2, berr(2,8)/2-berr(1,8)/2, berr(2,9)/2-berr(1,9)/2, berr(2,10)/2-berr(1,10)/2,0]];
    plot(E,I+30*(idx-1),'.-')
    plot(285:0.1:310, curve2(285:0.1:310)+30*(idx-1))
    end
    if Ts(idx)>66
        fo = fitoptions('Method','NonlinearLeastSquares',...
                'StartPoint',sp3);
        ft = fittype(eqn3,'options',fo);
        E=EsVHfit(idx,:);
    I=dataVHfit(idx,:);
    [curve2,gof2] = fit(E',I',ft);
    bsVH=[bsVH;[0,curve2.b1,curve2.b2,curve2.b3,0,0]];
    dsVH=[dsVH,curve2.c1];
    berr=confint(curve2,0.68);
    bsVHerr=[bsVHerr;[0, berr(2,4)/2-berr(1,4)/2, berr(2,5)/2-berr(1,5)/2, berr(2,6)/2-berr(1,6)/2,0,0]];
    plot(E,I+30*(idx-1),'.-')
    plot(285:0.1:310, curve2(285:0.1:310)+30*(idx-1))
    end
end

figure(40);
hold on
errorbar(Ts,bsVH-offsetsVH,bsVHerr,'diamond','Color',[0.93,0.69,0.13],'LineWidth',1,'MarkerFaceColor',[0.93,0.69,0.13],'MarkerSize',4)

ylim([290,306])
xlim([0,110])


figure(42);
hold on
plot(Ts,dsVH,'diamond','Color',[0.93,0.69,0.13],'LineWidth',1,'MarkerFaceColor',[0.93,0.69,0.13],'MarkerSize',4)
plot(Ts,dsHV,'o','Color',[0.85,0.33,0.10],'LineWidth',1,'MarkerFaceColor',[0.85,0.33,0.10],'MarkerSize',4)

xlim([0,110])

%% overplot


bs=[293,295,297,299,300.9,303]+0.5;
eqn='a1*exp(-(x-b1)^2/c1^2)+a2*exp(-(x-b2)^2/c1^2)+a3*exp(-(x-b3)^2/c1^2)+a4*exp(-(x-b4)^2/c1^2)+a5*exp(-(x-b5)^2/c1^2)+a6*exp(-(x-b6)^2/c1^2)+e*x+f';
sp=[10,30,50,30,10,10, bs(1), bs(2), bs(3), bs(4), bs(5), bs(6),.4,0,0];
bs=[293,295,297,299,300.9,303]+.54;
eqn2='a1*exp(-(x-b1)^2/c1^2)+a2*exp(-(x-b2)^2/c1^2)+a3*exp(-(x-b3)^2/c1^2)+a4*exp(-(x-b4)^2/c1^2)+a5*exp(-(x-b5)^2/c1^2)+e*x+f';
sp2=[10,30,50,30,10, bs(1), bs(2), bs(3), bs(4), bs(5),.4,0,0];
eqn3='a1*exp(-(x-b1)^2/c1^2)+a2*exp(-(x-b2)^2/c1^2)+a3*exp(-(x-b3)^2/c1^2)+e*x+f';
sp3=[10,20,20, bs(2)+.33, bs(3)+.3, bs(4)+.35,.4,0,0];
eqn4='a1*exp(-(x-b1)^2/c1^2)+e*x+f';
sp4=[20,298,5,0,0]
fo = fitoptions('Method','NonlinearLeastSquares',...
            'StartPoint',sp);
ft = fittype(eqn,'options',fo);
bsHV=[]
bsHVerr=[]
dsHV=[]
figure
hold on
for idx= 1:(length(Ts)-1):length(Ts)
    fo = fitoptions('Method','NonlinearLeastSquares',...
                'StartPoint',sp);
        ft = fittype(eqn,'options',fo);
    if Ts(idx)<45
        E=EsHVfit(idx,:);
        I=dataHVfit(idx,:);
        [curve2,gof2] = fit(E',I',ft);
        bsHV=[bsHV;[curve2.b1,curve2.b2,curve2.b3,curve2.b4,curve2.b5,curve2.b6]];
        dsHV=[dsHV,curve2.c1];
        berr=confint(curve2,0.68);
        bsHVerr=[bsHVerr;[berr(2,7)/2-berr(1,7)/2, berr(2,8)/2-berr(1,8)/2, berr(2,9)/2-berr(1,9)/2, berr(2,10)/2-berr(1,10)/2, berr(2,11)/2-berr(1,11)/2, berr(2,12)/2-berr(1,12)/2]];
        plot(E,I+30*(idx-1)*0-E*curve2.e-curve2.f,'.')
        plot(280:0.1:310, curve2(280:0.1:310)-(280:0.1:310)'*curve2.e-curve2.f+30*(idx-1)*0)
    end
    if Ts(idx)>66
        fo = fitoptions('Method','NonlinearLeastSquares',...
                'StartPoint',sp4);
        ft = fittype(eqn4,'options',fo);
        E=EsHVfit(idx,:);
    I=dataHVfit(idx,:);
    [curve2,gof2] = fit(E',I',ft);
    plot(E,(I-E*curve2.e-curve2.f),'o')
    plot(280:0.1:310, (curve2(280:0.1:310)-(280:0.1:310)'*curve2.e-curve2.f))
    end
end

%figure(44);
%errorbar(Ts,bsHV-offsetsHV,bsHVerr,'o','Color',[0.85,0.33,0.10],'LineWidth',1,'MarkerFaceColor',[0.85,0.33,0.10],'MarkerSize',4)
%ylim([290,310])
















